function [B,G] = get_navstk_bc(V,T,TE,ET,d,dof_map,cr,bdr,Re,caseNum)
%function [B,G] = get_navstk_bc(V,T,TE,ET,d,dof_map,cr,bdr,Re,caseNum)
% this is the boundary condition
nbdr = length(bdr);
Indx2 = zeros(nbdr*(2*d+1),1);
G = zeros(nbdr*(2*d+1),1);
row = 1;
G0 = 0;
% begin to treat boundary condition one by one bound edge
for k = 1:nbdr
    eg = bdr(k);
    tri = ET(eg,1);
    v1_loc = find(TE(tri,:)==eg);
    v2_loc = mod(v1_loc,3)+1;
    v3_loc = mod(v2_loc,3)+1;
    V1 = V(T(tri,v1_loc),:); 
    V2 = V(T(tri,v2_loc),:);  
    V3 = V(T(tri,v3_loc),:);  
    row_idx = row:row + 2*d;
    if v1_loc==2;
        v1_loc = -2;
    end
    line0 = cr_indices(0,d,v1_loc,cr);
    line1 = cr_indices(1,d,v1_loc,cr);
    [c1,c2] = navstk_first_layer(G0,V2,V3,V1,d,Re,caseNum); 
    G(row_idx) = [c1;c2];  
    Indx2(row_idx) = dof_map([line0;line1],tri);
    row = row + 2*d+1;
    G0 = c1(d+1);
end
N = length(Indx2);
ONE = ones(N,1);
Indx1 = (1:N)';
dim = max(max(dof_map))-min(min(dof_map))+1;
B = sparse(Indx1,Indx2,ONE,N,dim);
end

function [c1,c2] = navstk_first_layer(G0,V1,V2,V3,d,Re,caseNum)
d = d-1;
I = (0:d)';
J = d - I;
X = (J*V1(1) + I*V2(1))/d;
Y = (J*V1(2) + I*V2(2))/d;
[G,G1,G2] = boundary_function(X,Y,Re,caseNum); 
b3 = -G2*(V3(1)-V1(1)) + G1*(V3(2)-V1(2));
b2 = -G2*(V2(1)-V1(1)) + G1*(V2(2)-V1(2));
m = (d+1);
IM = diag(I)*ones(m,m);
JM = diag(J)*ones(m,m);
Mat = (IM/d).^(IM').*(JM/d).^(JM');
IF = gamma(I+1);
JF = gamma(J+1);
A = factorial(d)*ones(m,m)*diag(1./(IF.*JF));
Mat = A.*Mat;
D2 = 1/(d+1)*(Mat\b2);
B = diag(ones(d+1,1)) - diag(ones(d,1),-1);
D2(1) = D2(1) + G0;
c1 = B\D2;
c1 = [G0;c1];
c2 = 1/(d+1)*(Mat\b3) + c1(1:(d+1));

end

function line_dof = cr_indices(r,d,idx,pattern)
% to get the r's line of dofs of degree d parallel to edge i
% it is specially useful when treating smoothness condition and boundary
% condition! 
% the output parameter mat is the matrix P in my notes.
switch idx
    case 1
        line_dof = (((d+1-r)*(d-r)/2+1):((d+2-r)*(d+1-r)/2))';
    case -1
        line_dof = (((d+2-r)*(d+1-r)/2):-1:((d+1-r)*(d-r)/2+1))';
    case 2
        line_dof = pattern(r+1,1:d-r+1)';
    case -2
        line_dof = pattern(r+1,d-r+1:-1:1)';
    case 3
        line_dof = pattern(1:d-r+1,r+1);
    case -3
        line_dof = pattern(d-r+1:-1:1,r+1);
    otherwise
        line_dof = [];
end
%then put the corresponding value to a matrix;
% i = (1:d-r+1)';j = line_dof; s = ones(d-r+1,1);
% mat = sparse(i,j,s,d-r+1,(d+1)*(d+2)/2);
end

